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A direct TV-body model of core-collapse and core 
oscillations 
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ABSTRACT 

We report on the results of a direct TV-body simulation of a star cluster that started 
with N — 200 000, comprising 195 000 single stars and 5 000 primordial binaries. The 
code used for the simulation includes stellar evolution, binary evolution, an external 
tidal field and the effects of two-body relaxation. The model cluster is evolved to 
12 Gyr, losing more than 80% of its stars in the process. It reaches the end of the 
main core-collapse phase at 10.5 Gyr and experiences core oscillations from that point 
onwards - direct numerical confirmation of this phenomenon. However, we find that 
after a further 1 Gyr the core oscillations are halted by the ejection of a massive 
binary comprised of two black holes from the core, producing a core that shows no 
signature of the prior core-collapse. We also show that the results of previous studies 
with N ranging from 500 to 100 000 scale well to this new model with larger N. In 
particular, the timescale to core-collapse (in units of the relaxation timescale), mass 
segregation, velocity dispersion, and the energies of the binary population all show 
similar behaviour at different N. 

Key words: stars: evolution — globular clusters: general — galaxies: star clusters: 
general — methods: numerical — binaries: close — stars: kinematics and dynamics 



1 INTRODUCTION 

The increasing ability of the direct A^-body method to pro- 
vide reliable models of the dynamical evolution of star clus- 
ters has closely mirrored increases in computing power (Heg- 
gie 2011). The community has progressed from small- N 
models performed on workstations (e.g. von Hoerner 1963; 
McMillan, Hut & Makino 1990; Giersz & Heggie 1997) to 
models of old open clusters and into the N ~ 100 000 regime 
(Baumgardt & Makino 2003) by making use of special- 
purpose GRAPE hardware (Makino 2002). Software ad- 
vances over a similar timeframe have produced sophisticated 
codes (Aarseth 1999; Portegies Zwart et al. 2001) that in- 
crease the realism of the models by incorporating stellar 
and binary evolution, binary formation, three-body effects 
and external potentials. As a result, Af-body models have 
been used in numerous ways to understand the evolution 
of globular clusters (GCs: Vesperini & Heggie 1997; Baum- 
gardt & Makino 2003; Zonoozi et al. 2011), even though the 
best models still only touch the lower end of the GC mass- 
function (see Aarseth 2003 and Heggie & Hut 2003 for a 
more detailed review of previous work). At the other end 
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of the spectrum, Monte Carlo (MC) models have proven ef- 
fective at producing dynamical models of N > 10 6 particles 
(Giersz & Heggie 2011). These models have shown that clus- 
ters previously defined as non-core-collapse can actually be 
in a fluctuating post-core-collapse phase (Heggie & Giersz 
2008). In practice the two methods are complimentary with 
MC informing the more laborious A-body approach (such as 
refining initial conditions) and A-body calibrating aspects 
of MC. 

In this paper we present an A^-body simulation of star 
cluster evolution that begins with N = 200 000 stars and 
binaries. This extends the A^ parameter space covered by 
direct A-body models and performs two important func- 
tions. Firstly it provides a new calibration point for the MC 
method - this statistical method is increasingly valid for in- 
creasing N so calibrations at higher N are more reliable. 
It also allows us to further develop our theoretical under- 
standing of star cluster evolution and investigate how well 
inferences drawn from models of smaller N scale to larger 
values. The latter is the focus of this current paper. A good 
example of the small-A models that we wish to compare 
with is the comprehensive study of star cluster evolution 
presented by Giersz & Heggie (1997) using models that in- 
cluded a mass function, stellar evolution and the tidal field 
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of a point-mass galaxy, albeit starting with 500 stars in- 
stead of 200 000. More recent examples for comparison in- 
clude Baumgardt & Makino (2003) and Kupper et al. (2008). 
We were also motivated to produce a model that exhibited 
core-collapse close to a Hubble time without dissolving by 
that time. What we find when interpreting this model is 
that much of the behaviour reported previously for smaller 
JV-body models stands up well in comparison but that the 
actions of a binary comprised of two black holes (BHs) pro- 
vides a late twist to the evolution of the cluster core. 

In Section 2 we describe the setup of the model. This is 
followed by a presentation of the results in Sections 3 to 7 
focussing on general evolution (cluster mass and structure), 
the impact of the BH-BH binary, mass segregation, velocity 
distributions and binaries (binary fraction and binding en- 
ergies). Throughout these sections the results are discussed 
and compared to previous work where applicable. Then in 
Section 8 we specifically look at how the evolution timescale 
of the new model compares to findings presented in the past. 



2 THE MODEL 

For our simulation we used the NB0DY4 code (Aarseth 1999) 
on a GRAPE-6 board (Makino 2002) located at the Ameri- 
can Museum of Natural History. NB0DY4 uses the 4th-order 
Hermite integration scheme and an individual timestep al- 
gorithm to follow the orbits of cluster members and invokes 
regularization schemes to deal with the internal evolution 
of small- JV subsystems (see Aarseth 2003 for details). Stel- 
lar and binary evolution of the cluster stars are performed 
in concert with the dynamical integration as described in 
Hurley et al. (2001). 

The simulation started with 195 000 single stars and 

5 000 binaries. We will refer to this as the K200 model. The 
binary fraction of 0.025 is guided by the findings of Davis et 
al. (2008) which indicated a present day binary fraction of 
~ 0.02 for the globular cluster NGC 6397, measured near the 
half-light radius of the cluster. As shown in Hurley, Aarseth 

6 Shara (2007) and discussed in Hurley et al. (2008), this 
can be taken as representative of the initial binary fraction 
of the cluster. Thus we adopted this value for our model. 
Validation of the binary fraction approach will be provided 
in Section 7. 

Masses for the single stars were drawn from the initial 
mass function (IMF) of Kroupa, Tout & Gilmore (1993) be- 
tween the mass limits of 0.1 and 50Mq. Each binary mass 
was chosen from the IMF of Kroupa, Tout & Gilmore (1991), 
as this had not been corrected for the effect of binaries, and 
the component masses were set by choosing a mass-ratio 
from a uniform distribution. In NB0DY4 we assume that all 
stars are on the zero-age main sequence when the simula- 
tion begins and that any residual gas from the star formation 
process has been removed. A metallicity of Z — 0.001 was 
set for all stars. The orbital separations of the 5 000 primor- 
dial binaries were drawn from the log-normal distribution 
suggested by Eggleton, Fitchett & Tout (1989) with a peak 
at 30 au and a maximum of 100 au. Orbital eccentricities of 
the primordial binaries were assumed to follow a thermal 
distribution (Heggie 1975). 

For the tidal field of the parent galaxy we have used 
the point-mass galaxy approach with the model cluster on 



a circular orbit at R gc — 3.9 kpc with an orbital velocity 
of V g = 220 km s" 1 . In setting R gc we have been primarily 
guided by a desire for the cluster to have its moment of core- 
collapse between 10 — 13 Gyr. Previous experience suggested 
that R gc ~ 4 kpc would provide this for a model starting 
with JV = 200 000. 

We used a Plummer density profile (Plummer 1911; 
Aarseth, Henon & Wielen 1974) and assumed the stars and 
binaries are in virial equilibrium when assigning the initial 
positions and velocities. The Plummer profile formally ex- 
tends to infinite radius so in practice a cut-off at a radius of 
~ 10 rh is applied, where is the half- mass radius. This is 
to avoid rare cases of large distance in the initial distribu- 
tion. The tidal field sets a tidal radius according to: 




where G is the gravitational constant and M is the clus- 
ter mass (see Giersz & Heggie 1997). We chose the JV-body 
length-scale of our model so that the outermost star of the 
initial model sits at r t . This reflects the expansion expected 
when gas leftover from star formation is removed from the 
potential well (which we do not model), hence our initial 
model should be more radially extended than a compact 
protocluster. Stars were removed from the simulation when 
their distance from the density centre exceeds twice that of 
the tidal radius of the cluster. 

With these choices the initial parameters of the K200 
model were M = 100 067.4 Mq , = 4.7pcandr t = 35.8 pc. 
The half-mass relaxation timescale of the initial model was 
t A ~ 1 115 Myr. 

Much of the behaviour of this model will be compared 
to that reported by Hurley et al. (2008) for a model that 
started with 95 000 single stars and 5 000 binaries. This will 
be referred to as the K100 model. The K100 model was 
placed on a circular orbit about a point-mass galaxy at a 
radial distance of i? gc = 8.5 kpc. It had the same number of 
binaries as the K200 model but twice the binary fraction. 
The parameters of the stars and binaries were set up in the 
same manner as described above for the K200 model. 



3 GENERAL EVOLUTION 

In Figure Q] we look at the evolution of the total cluster 
mass with time for the K200 model. This simulation was 
stopped for analysis at 12 Gyr having satisfied the goal of 
providing a post-core-collapse model (see below) at the ap- 
proximate age of a Milky Way globular cluster. At this point 
the model cluster has lost 83% of its initial mass. In terms of 
stars remaining at 12 Gyr the K200 model has JV = 31 713 
comprised of 30 835 single stars and 878 binaries. 

We include the evolution of the K100 model from Hur- 
ley et al. (2008) in Figure [T] and see that at 12 Gyr the two 
model clusters have the same amount of mass remaining (af- 
ter starting with a factor of two difference). For comparison 
we also show in Figure [T] a model that started with 100 000 
stars on the same orbit as the K200 model. As expected this 
model does not last for a Hubble time and is completely dis- 
solved after about 9 Gyr. The slope in the mass-age plane 
is similar to that of the K200 model and distinct from that 
of the K100 model with J? gc = 8.5 kpc. An investigation of 
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Figure 1. Evolution of the total cluster mass as a function of 
time for the N = 200 000 model orbiting at R gc = 3.9 kpc (solid 
line), the Hurley et al. (2008) N = 100 000 (KlOO) model orbiting 
at Rg C = 8.5 kpc (dashed line) and a N = 100 000 model orbiting 
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Figure 2. Evolution of the A-body core radius (+ symbols), the 
half-mass radius (lower solid line) and tidal radius (upper solid 
line) for the N = 200 000 model. All radii are three-dimensional. 
The vertical dotted line at 10.5 Gyr denotes the time we have 
identified with the end of the initial core-collapse phase. 



the mass-loss rates and dissolution times of star clusters as 
a function of orbit within the Galaxy will be the subject of 
another paper (Madrid et al. 2012). 

Figure [2] shows the behaviour of the core radius, the 
half-mass radius and the tidal radius as the K200 model 
evolves. Both the half-mass and core radii show an initial 
increase corresponding to stellar evolution mass-loss from 
massive stars, which are mostly found in the inner regions 
of the cluster. The half-mass radius then plateaus before ap- 
pearing to follow the decreasing trend of the tidal radius at 
later times. At 12 Gyr we have ~ 3pc which is compara- 
ble to the average effective radius found for globular clusters 
(e.g. Jordan et al. 2005). 

The core-radius shows a deep minimum at 10.5 Gyr 
which we identify as the moment that the initial core- 
collapse phase ends. This is based purely on inspection of 
Figure [2] noting that the same method - looking for the 
first deep minimum of the density- or mass-dependent cen- 
tral radius - has been commonly employed in the past (e.g. 
Baumgardt & Makino 2003; Hurley et al. 2004; Kiipper et al. 
2008). At this point the core density is 5 200 stars pc -3 , in- 
creased from an initial value of 700 stars pc -3 . Subsequently 
the core fluctuates markedly, corresponding to post-core- 
collapse oscillations highlighted by Heggie & Giersz (2009). 
However, we note that the core exhibits fluctuating be- 
haviour leading up to the moment of core-collapse as well. 

The core radius in Figure [2] is the density radius com- 
monly used in A^-body simulations (Casertano & Hut 1985), 
calculated from the density weighted average of the distance 
of each star from the density centre (Aarseth 2003). Fol- 
lowing Heggie & Giersz (2009) we have also looked at the 
dynamical core radius used in their Monte Carlo models, 
calculated as r^ D = 3wfo/ (47rGp2o) where V20 is the mass- 



weighted central velocity dispersion and P20 is the central 
density, both calculated from the innermost 20 stars. We 
find that r c and r Ci D cover a similar range at all times. In 
particular, for the period 10.5 — 11.5 Gyr, i.e. 1 Gyr after 
core- collapse, both oscillate between 0.1 to 0.6 pc for the 
majority of the time. For reference, the radius containing 
the inner 1% of the cluster mass is ~ 0.3 pc over this pe- 
riod and is thus a good proxy for the average core radius. 
Also following Heggie & Giersz (2009) we have used the 
autocorrelation method to determine if there is any clear 
periodicity in the fluctuations of the A^-body core radius in 
the 1 Gyr subsequent to core-collapse. This showed a pe- 
riod of about 120 Myr: greater than the crossing time (few 
Myr) and less than the relaxation time (300 — 400 Myr). In 
comparison, Heggie & Giersz (2009) reported an oscillation 
period of ~ 400 Myr for their A^-body modeQ, although this 
contained more A" than our post-collapse K200 model and 
correspondingly had a higher half-mass relaxation timescale 
of ~ 700 Myr. 

It is of interest to examine how previous A-body results 
reported for smaller N models hold up in comparison to our 
new model. The KlOO model of Hurley et al. (2008) reached 
a similar core radius at the end of core-collapse as did our 
K200 model (r c ~ 0.1 pc), albeit at a later time (16 Gyr com- 
pared to 10.5 Gyr). As noted above, the average value of r c 
near core-collapse is similar to the 1% Lagrangian radius for 
the K200 model. For the KlOO model r c evolves similarly 
to the 2% Lagrangian radius (although it does dip down to 
the 1% radius on occasion) and for the A = 500 models of 
Giersz & Heggie (1997) r c is at the 5% Lagrangian radius 

1 This A-body model was evolved for 1 Gyr starting from a post- 
collapse Monte Carlo model at an age of 12 Gyr. 
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or greater. Thus it seems that with increasing TV the depth 
of core-collapse increases relative to the cluster mass distri- 
bution. If we look instead at scaled quantities, particularly 
the evolution of the ratio r c /rh as a function of age scaled 
by the half-mass relaxation timescale, we find that the K200 
and K100 simulations track each other very well. The ratio 
starts at ~ 0.4 and steadily decreases to an average value 
of ~ 0.1 at the end of core-collapse. Giersz & Heggie (1997) 
find r c /rh ~ 0.15 for their models at core-collapse, noting 
that they use r Cj o rather than the Casertano & Hut (1985) 
definition and that latter is about twice as large in the post- 
collapse phase for their models (we found that the two gave 
similar average values). McMillan (1993) performed models 
with TV = 1 024, primordial binaries and a tidal field. For 
these models r c jr^ stabilized at about 0.1 in good agree- 
ment with prior models of isolated clusters (see McMillan 
1993 for details). It therefore seems that this can be taken 
as a reliable value for clusters at the end of core-collapse 
(although see the next section for mention of some unusual 
cases) . 



4 THE LATE ACTION OF A BLACK-HOLE 
BINARY 

A remarkable feature in Figure is the sharp change in the 
behaviour of the core at ~ 11.5 Gyr, the last deep minimum, 
when the size of the core suddenly increases and evolves 
steadily from that point onwards. This change is related to 
an interaction within the core involving a binary comprised 
of two BHs. 

The binary in question is non-primordial. Each BH 
formed from a massive single star within the first 10 Myr 
of evolution with masses of 23.6 and 17.3 Mq, respectively. 
The two BHs formed a binary at 1 800 Myr in a four-body 
interaction, initially with a very high eccentricity and long 
orbital period of 10 5 d. It resided in the core for the majority 
of its lifetime and suffered a series of perturbations and in- 
teractions that saw the eccentricity vary between 0.2 to 0.95 
and the orbital period reduced to 13 d at 11 550 Myr. At that 
time the BH-binary becomes embroiled in a strong interac- 
tion with a binary comprised of two main-sequence stars 
(masses of 0.8 and 0.6 Mq). The second binary is broken-up 
and the two main-sequence stars are ejected rapidly from 
the cluster (velocities of 96 and 154 kms -1 ). This causes a 
recoil of the BH-binary which leaves the core and then the 
cluster entirely (10 Myr later with a velocity of 5kms _1 ). 
The domination of the central region of the cluster by this 
BH-binary and its subsequent ejection are similar to the 
processes described by Aarseth (2012). 

The sudden loss of mass from the core - the average 
mass drops by 30% (see next section) - combined with the 
rapid ejection of the two main-sequence stars causes the core 
to expand. We see that after this event the core radius does 
start to decrease once more but without fluctuations. Thus, 
the influence of one strong interaction involving a massive 
binary has halted the core oscillation process. Compared to 
the point that we identified as the end of the initial core- 
collapse phase the core radius has increased by a factor of 
about six. The structure of globular clusters is often quan- 
tified by the concentration parameter c = log r t /r c (King 
1966). Milky Way GCs exhibit a range of c values (Harris 
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Figure 3. Average mass in different Lagrangian regions (0-1%, 
1-10%, 10-50%, 50-100%) for the N = 200 000 model. The vertical 
dotted line marks the end of the initial core-collapse phase as in 
Figure [2] 

1996) with the most obvious core collapse examples having 
c = 2.5 but with c > 2 generally taken as indicative of a 
high-density cluster or a possible core-collapse cluster (Ma- 
teo 1987). At the end of the core-collapse phase our K200 
model has c ~ 2.3 and this decreases to c ~ 1.5 after the BH- 
binary is ejected from the core. Thus, the cluster would not 
be expected to appear as a core-collapse cluster if observed 
at this point. 

Hurley (2007) showed that the presence of a long-lived 
BH-BH binary in the core, with both BHs being of stellar 
mass, could significantly increase the r c /rh ratio of a model 
with 100 000 stars. The BH-BH binary in our TV = 200 000 
model has not produced a similar inflation of the ratio. 
Mackey et al. (2007) performed TV-body simulations with 
TV ~ 100 000 in which they retained ~ 200 stellar mass BHs. 
They found that the BHs formed a dense core in which in- 
teractions were common and BHs could be ejected from the 
cluster, leading to a significantly inflated core radius. It is 
our intention in the near future to look at a wide range of TV- 
body simulations and document in detail the statistics and 
outcomes of BH-BH binaries in the cores of model clusters. 
This will include fitting King (1966) models to the density 
profiles of the model clusters so as to properly calculate the 
concentration parameter rather than using TV-body values 
as we have done here in our preliminary analysis. 



5 MASS SEGREGATION 

Figure [3] looks at how the average stellar mass behaves for 
the K200 model, focussing on four different Lagrangian re- 
gions: a central volume that encompasses the inner 1% of the 
cluster mass, a central shell that lies between radii enclosing 
1% and 10% of the cluster mass, an intermediate shell that 
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Figure 4. As for Figure [3] but now showing the velocity disper- 
sion. 



lies between the 10% and 50% Lagrangian radii, and an outer 
shell that includes all stars beyond the 50% Lagrangian ra- 
dius. Note that binaries are included and are assumed to 
be unresolved. The average stellar mass throughout the en- 
tire cluster is O.5M0 at the start of the simulation - there 
is no primordial mass segregation - and drops initially in 
all regions owing to stellar evolution mass-loss of massive 
stars for the first ~ 100 Myr. We then see that the effect of 
mass segregation driven by two-body encounters takes over, 
causing an increase in the average stellar mass in the inner 
regions and a corresponding decrease in the outer region. By 
the time that one half-mass relaxation timescale has elapsed 
(~ 1 000 Myr) there is a clear distinction between the aver- 
age mass in each of the regions. In the central regions the 
average stellar mass continues to increase up to the end of 
core-collapse and then flattens post-collapse (until the ejec- 
tion of the BH-BH binary). The value in the very centre is 
noisy owing to a smaller number of objects and the cycling 
of these objects in and out of the region. We see a marked 
increase in this central value as the cluster gets closer to 
the end of core-collapse (note the correlation with r c in Fig- 
ure [2} . The decreasing average mass in the outer region is 
gradually arrested by the effect of the external tide which 
preferentially removes the lower-mass stars that have been 
pushed out to this region. We see that from ~ 5 000 Myr on- 
wards the average mass in the outer region is now increasing 
and that this effect is also felt in the intermediate region. 

If we instead look at the behaviour in two-dimensional 
projected regions we find that the average mass is the same 
in the two outermost regions but drops by about 5% in the 
1-10% region and 10-15% in the central region (after the 
first Gyr). 

Giersz & Heggie (1997) look at the evolution of aver- 
age mass in different regions, as we have done in Figure [3] 
They see very similar behaviour which can be summarised 
as: (i) a sharp increase of the average mass within the 1% 
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Figure 5. As for Figure [3] but now showing the anisotropy pa- 
rameter. 

Lagrangian region towards core- collapse; (ii) a decrease in 
the outer regions that flattens out with time and then in- 
creases at late times; and (iii) similar trends in the interme- 
diate regions. However, in their N = 500 models they find 
that the timescale for the initial phase of mass segregation 
is about the same as the core-collapse timescale, whereas 
we find that mass segregation is fully established well be- 
fore core-collapse. The iV = 131 072 models of Baumgardt 
& Makino (2003) agree with our K200 model in that respect 
and with the general behaviour, including the flattening of 
the average mass in the central regions post-collapse. 



6 VELOCITY DISTRIBUTIONS 

In Figure [4] we show the velocity dispersion for the same 
Lagrangian regions as in Figure [3] All regions show a rapid 
initial decrease owing to an overall expansion of the cluster. 
This is followed by a more gradual decrease as the cluster 
evolves towards core-collapse, with all regions declining in 
an almost homologous manner. As the model cluster nears 
the end of core-collapse there is a pronounced upturn in the 
velocity dispersion within the 1% radius. This is also seen 
out to the 10% radius and even at the half-mass radius, al- 
though to a much smaller extent. Note that the behaviour 
for two-dimensional velocities in projected Lagrangian re- 
gions is the same but with values that are typically 20% 
less. 

The main features of our new model mirror those in the 
smaller- N models of Giersz & Heggie (1997). These features 
are the following: the values within the 1% and 10% regions 
overlap; the velocity dispersion in the outer regions is clearly 
lower than for these inner regions; and the values for the 50% 
region are much closer to those of the inner regions than the 
outer regions. 

Another aspect of velocity to consider is the anisotropy. 
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We have chosen to define this as the ratio of the mean square 
transverse to radial velocity components, < v% > / < > 
(as in Giersz & Heggie 1997) and the result is shown in Fig- 
ure [5] In three dimensions this is equal to two for isotropy, 
greater than two for a tangentially anisotropic distribu- 
tion and less than two for a radially anisotropic distri- 
bution. An alternative is to use the anisotropy parameter 
P = 1- < vt > /2 < Vr > (Binney & Tremaine 1987; Baum- 
gardt & Makino 2003; Wilkinson et al. 2003) where /3 = 
is isotropic, f3 < is tangentially anisotropic and f3 > is 
radially anisotropic. We see from Figure[5]that the inner re- 
gions of the cluster remain close to isotropic throughout the 
evolution while the outer region develops an increasing tan- 
gential anisotropy over the first ~ 7Gyr and then flattens 
out to a constant value for the remainder of the evolution. 
This overall behaviour is similar to that observed by Baum- 
gardt & Makino (2003) in their models. It is also similar to 
the behaviour for the K100 model, although the degree of 
anisotropy in the outer region is about 30% less than in the 
K200 model. Also, because the K100 model is evolved well 
past core-collapse it is possible to see that the anisotropy is 
reduced post-collapse and tends towards isotropy in the final 
stages of evolution (as was noted by Baumgardt & Makino 
2003). The tangential anisotropy in the outer region is con- 
trary to previous findings of radial anisotropy for smaller- iV 
models (Giersz & Heggie 1997; Wilkinson et al. 2003). Tan- 
gential anisotropy can be explained by stars expelled from 
the central regions on radial orbits preferentially escaping 
from the cluster at the expense of stars on tangential orbits 
which find it harder to escape. In the very central region 
(within the 1% Lagrangian radius) the anisotropy parame- 
ter is very noisy and fluctuates between radial and tangential 
anisotropy. However, the average behaviour is isotropy. 



7 BINARIES 

Hurley, Aarseth & Shara (2007) documented the evolution 
of binary fraction with time for a range of TV-body models 
covering N = 24 000 to 100 000 and primordial binary frac- 
tions from 0.05 to 0.5. In all cases they found that the binary 
fraction within the cluster core and the 10% Lagrangian ra- 
dius increases as the cluster evolves towards core-collapse 
while the overall binary fraction of the cluster stays close 
to the primordial value. Figure |S] shows the same evolution 
for the K200 model. Similar results are seen in that the core 
binary fraction increases markedly as the cluster evolves, 
particularly towards core-collapse when it becomes a fac- 
tor of six or more greater than the primordial value, while 
the overall binary fraction does stay close to the primordial 
value (although the inclusion of a large proportion of ini- 
tially wide binaries would lead to an initial decrease). This 
gives added confidence in the results of Hurley, Aarseth & 
Shara (2007) although models of greater N (and density) 
are still required before the situation for the majority of the 
Milky Way globular clusters can be firmly established. 

Davis et al. (2008) measured the binary fraction 
amongst main-sequence stars near the half-mass radius of 
NGC6397 and found it to be a few per cent. This was 
taken as representative of the primordial binary fraction of 
NGC 6397 by leaning on the result from Hurley et al. (2008) 
showing that the binary fraction measured near the half- 
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Figure 6. Evolution of the binary fraction for the N = 200 000 
model. Shown are the binary fraction of the entire cluster (dashed 
line), within the 10% Lagrangian radius (black dotted line) and 
within the core (red solid line). Also shown is the binary fraction 
of main-sequence stars and binaries near the half-mass radius 
(red dotted line: which closely follows the binary fraction of the 
entire cluster) . Compare with Figure 3 of Hurley, Aarseth & Shara 
(2007). The vertical dotted line once again marks the end of the 
initial core-collapse phase as in Figure [2] 



mass radius of a cluster can be taken as a good indication 
of the primordial binary fraction. This is reinforced in Fig- 
ure [6] for our K200 model where we show that the binary 
fraction of main-sequence stars and binaries near the half- 
mass radius of the cluster stays at roughly the same value 
throughout the evolution. 

Next we look at the energy in binaries for the K200 sim- 
ulation in Figure [7]and for the K100 simulation in Figure [5] 
Here we see that for the 10— 100% Lagrangian mass regions 
the values and behaviour for the two simulations are very 
similar. As expected the average energy per binary increases 
as we move inwards towards the cluster centre, although the 
values within the 1-10% and 10-50% regions are converging 
towards the end of the simulation. We note that the binary 
binding energy is calculated in arbitrary yet physical units 
of Mq / Rq to enable direct comparison between the internal 
energies of the two binary populations. Conversion to cluster 
units of fcT based on the average kinetic energy of the cluster 
stars (e.g. McMillan 1993) increases the K100 values by a 
factor of three compared to the K200 values, i.e. the binaries 
have relatively more internal energy in the simulation with 
less stars. Sharp dips in the binding energy occur through- 
out the simulations when binaries escape (as documented 
by Kiipper, Kroupa & Baumgardt 2008) although these are 
less evident in Figure [5] owing to less frequent sampling. 

The energy per binary in the inner 1% region by mass 
starts off similarly for both simulations, oscillating about 
the 1-10% values. However, the influence of the BH-binary 
in the K200 simulation from ~ 5 Gyr onwards is clear to 
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Figure 7. Average energy per binary in different Lagrangian 
regions (0-1%, 1-10%, 10-50%, 50-100%) for the N = 200 000 
model. The vertical dotted line denotes the time identified with 
the end of the initial core-collapse phase for this model. 




t/Myr 

Figure 8. As for Figure [7] but now for the K100 model. The 
vertical dotted line denotes the time identified with the end of 
the initial core-collapse phase for this model. 

see, increasing the average energy by almost two orders of 
magnitude before a sharp drop when the binary is removed 
from the core. 

Noticeable increases in binding energy are also evident 
at other times and arise from the creation of tight binaries 
by various means. For example, the spike in the 50-100% 
Lagrangian region at ~ 4 300 Myr in Figure [Jj results from 



Table 1. Comparison of core-collapse, i cc , and relaxation 
timescales, t^, for the K200 and K100 simulations. Note that 
the averaged tjj, is calculated to 11.5 Gyr and 18Gyr for the 
K200 and K100 simulations, respectively. Also shown is the 
time for the cluster to lose half of the initial mass, t 1 / 2 - 



quantity 


K200 


K100 


*i/ 2 /Myr 


4042 


5415 


t cc /Myr 


10500 


16000 


t r h/Myr (initial) 


1115 


1405 


t rh /Myr (*i/ 2 ) 


1480 


1910 


t rh /Myr (cc) 


415 


500 


trh/Myr (average) 


762 


1340 



a primordial binary in which common-envelope evolution 
creates a short-period binary comprised of two white dwarfs 
which subsequently come into contact and merge. Thus the 
binary dominates the energy in this region for a short time 
and then disappears. We see a more persistent increase in 
the average energy within the 1% region at ~ 4 500 Myr 
in Figure [7] This is caused when a main-sequence star and 
a black hole form a short-period binary via an exchange 
interaction. The binary survives for about 400 Myr in the 
cluster core before the main-sequence star is consumed by 
the black hole. 

It is worthwhile to ask if the signatures of these short- 
lived energetic binaries be observed? The binary that re- 
sulted in the merger of two WDs did not produce a WD 
with a mass in excess of the Chandrasekhar limit so could 
not be a potential type la supernova (see Shara & Hurley 
2002). However, the resulting WD will be relatively massive, 
hot and young. Thus it could be expected to remain one of 
the brightest WDs in the cluster for several Gyrs and would 
be easily observed. The other binary mentioned resulted in 
a main-sequence star being tidally disrupted and swallowed 
by a black- hole of ~ 20 Mq . This could be expected to give 
rise to a burst of X-rays and perhaps gamma-rays, along the 
lines of Burrows et al. (2011: albeit for a supermassive black 
hole). 



8 SCALING OF PREVIOUS TIMESCALE 
RESULTS 

Timescales for star cluster evolution are important to under- 
stand, with the time until core-collapse and the time until 
dissolution being quantities of interest (e.g. Gnedin & Os- 
triker 1997). Furthermore, the scaling of these with N or re- 
lated properties of clusters/models is necessary (Baumgardt 
2001), particularly while direct models of globular clusters 
remain out of reach. 

In Table [1] we summarize the significant timescales for 
the K200 and K100 simulations: the time for half of the 
initial mass to be lost, the time until the end of the core- 
collapse phase and the half-mass relaxation time at various 
key points. 

Kiipper et al. (2008) presented a range of open clus- 
ter models in a steady tidal field and found that the core- 
collapse time scaled by the initial half-mass relaxation time 
ranged from 17 for clusters starting well within their tidal 
radii (smaller t r h,o) to 9 for clusters that fill their tidal radii 
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from the start (larger t r h,o). Our K200 model and the K100 
model of Hurley et al. (2008) start by filling their tidal ra- 
dius and have i C c/£rh,o ~ 10 which is in good agreement with 
Kiipper et al. (2008). 

Baumgardt (2001) looked at the scaling of A-body 
models using simulations of up to 16 384 equal-mass stars. 
In this work the time for a cluster to lose half of its mass, 
ti/2, was taken as an indication of the cluster lifetime with 

3/4 

ti/2 * r h shown to be the appropriate scaling. Baumgardt 
& Makino (2003) subsequently used their larger-A models 
(up to 131 072 stars) to once again show that scaling the dis- 

3/4 

solution time as is appropriate, this time for multi-mass 
models that included stellar evolution. However, we see from 
Table Q] that the half-mass relaxation time varies consider- 
ably across the lifetime of a cluster and it is not immediately 
clear which value to use when scaling timescales. For an ob- 
served cluster it will be the value at the current age of the 
cluster. A fairer comparison would be to use the average half- 
mass relaxation time across the lifetime of the cluster, i r h,av 
- of course this is not known for an observed cluster. For the 
K200 and K100 models we find that tcc/t^ 4 is the same for 
both simulations, so the agreement for this key timescale is 
in excellent agreement with the previous suggestion. We also 

3/4 

find that i]/ 2 scales quite well with t r £ . However, if we in- 
stead use £^ 4 q then we do not find good agreement for either 
the core-collapse or half-life timescales. Thus our agreement 
with the scalings found in previous works is dependent on 
which t T h is used. 



9 SUMMARY 

We have presented an A-body model that started with 
200 000 stars and binaries, evolves to the moment of core- 
collapse at 10.5 Gyr and has ~ 30 000 stars remaining at 
12 Gyr. We have used our direct A-body model to confirm 
the post-core-collapse fluctuations described in the Monte 
Carlo model of Heggie & Giersz (2008) and the hybrid A- 
body/MC approach of Heggie & Giersz (2009). We have also 
shown that these fluctuations can be halted by the ejection 
of a dominant BH-binary from the core. This produces a core 
that shows no sign that it has previously evolved through 
core- collapse. 

We have looked at how the results of previous works 
compare to a model of larger A and find good agreement 
provided that appropriate scalings are used (such as the 
core-radius to half- mass radius ratio at core-collapse). In 
terms of raw values some variations exist: the core radius 
at core-collapse reaches deeper in to the mass distribution 
for larger A, for example. The behaviour of quantities such 
as average stellar mass and velocity dispersion have been 
documented and the general behaviour matches expecta- 
tions from earlier models. Looking at time scales such as 
the time to core-collapse and the dissolution time we also 
find agreement with scaling relations previously reported in 
the literature, however this is dependent on which value of 
the half-mass relaxation timescale is used. In particular, the 
scaling of dissolution time with tj?^ 4 reported by Baumgardt 
(2001) could be reproduced provided that the average £ r h 
was used and not the initial i r h. 

The A — 200 000 simulation reported in this work took 



the best part of a year on a GRAPE-6 board to complete. It 
continues the gradual increase of A used in realistic A-body 
models from the A = 500 model of Giersz & Heggie (1997) 
to the A = 131 072 model of Baumgardt & Makino (2003). 
However, with only ~ 20 000 Mq remaining in our model at 
an age of 11.5 Gyr we are still only touching the lower end 
of the globular cluster mass-function. This is after consid- 
erable effort. In particular, we are still some way from the 
goal of a full million-body model of a globular cluster (Heg- 
gie & Hut 2003). How can we push forward to reach that 
goal? The shift towards graphics processing units (GPUs) as 
the central computing engine for A-body codes, combined 
with sophisticated software development, offers hope (Nita- 
dori & Aarseth 2012). Simulations of 100 000 stars can be 
performed comfortably on a single-GPU (Hurley & Mackey 
2010; Zonoozi et al. 2011) and the introduction of multiple- 
GPU support will likely make simulations of the type pre- 
sented here commonplace in the near future. Further hard- 
ware advances and a revisiting of efforts to parallelize direct 
A-body codes (Spurzem 1999) will also aid the push towards 
greater A. 

In a follow-up paper we will conduct a full investigation 
of the stellar and binary populations of our model. It is also 
our intention to make model snapshots - saved at frequent 
intervals across the lifetime of the simulation - available for 
others to observe and analyse. These can be obtained by 
contacting the authors. At the beginning of this paper we 
indicated that an important function of a large-A model 
would be to aid in the calibration of the Monte Carlo tech- 
nique. This is currently underway (Giersz et al. 2012) and 
will include a direct comparison of A-body and MC models 
starting from the same initial conditions. 
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